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Homogeneous nucleation of dislocations as bifurcations in a peri- 
odized discrete elasticity model 
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Abstract. - A novel analysis of homogeneous nucleation of dislocations in sheared two- 
dimensional crystals described by periodized discrete elasticity models is presented. When the 
crystal is sheared beyond a critical strain F — F c , the strained dislocation- free state becomes 
unstable via a subcritical pitchfork bifurcation. Selecting a fixed final applied strain Ff > F c , 
different simultaneously stable stationary configurations containing two or four edge dislocations 
may be reached by setting F = Fft/t r during different time intervals t r . At a characteristic time 
after t r , one or two dipoles are nucleated, split, and the resulting two edge dislocations move in 
opposite directions to the sample boundary. Numerical continuation shows how configurations 
with different numbers of edge dislocation pairs emerge as bifurcations from the dislocation-free 
state. 



Introduction. — Homogeneous nucleation of disloca- 
tions is observed in different processes such as nanoinden- 
tation experiments [1,2], heteroepitaxial crystal growth 
[3,4], indentation experiments in colloidal crystals [5] or 
soap bubble raft models [6]. Homogeneous nucleation of 
dislocations occurs in a perfect crystal and is therefore ex- 
pected to have a much higher activation energy than het- 
erogeneous nucleation at defect sites such as step edges. 
Different types of calculations have been used to interpret 
homogeneous nucleation of dislocations in different situ- 
ations, ranging from atomistic simulations to continuum 
mechanics interpretations or combinations thereof [7]. In 
all reliable nucleation criterion is needed to cap- 

ture the nature of nucleated defects and the time and place 
at which such defects appear. 

In this letter, we tackle homogeneous nucleation of dis- 
locations as a bifurcation problem in discrete elasticity. 
Periodized discrete elasticity models of dislocations regu- 
larize in a natural way the singularities of the stress field 
at the dislocation lines that appear in the theory of elas- 
ticity, and they also allow the sliding motion of crystal 
planes typical of dislocation gliding [8,9]. Simplified ver- 
sions of these models have been used to analyze disloca- 
tion depinning and motion at the Peierls stress in a precise 



manner [10]. Here we show that simple periodized discrete 
elasticity models are able to describe homogeneous nucle- 
ation of dislocations by shearing an initially undisturbed 
dislocation-free lattice. While molecular dynamics pro- 
duces nucleation of dislocations and the Peierls stress, it is 
considerably costlier than periodized discrete elasticity. In 
more efficient computer codes of discrete line dislocation 
dynamics, nucleation of dislocations and the Peierls stress 
have to be added to the code [7]. Moreover, simple mod- 
els allow for a more detailed analysis and interpretation of 
results than either these computer- intensive methods. 

The idea behind periodized discrete elasticity models 
is to discretize space according to the nodes of the crys- 
tal lattice. Then the potential energy of the crystal is 
constructed as follows. The usual strain energy density 
of linear elasticity is obtained from the tensor of elastic 
constants and the strain tensor which is the symmetrized 
gradient of the displacement vector. In periodized dis- 
crete elasticity models of simple cubic crystals, the gra- 
dient of the displacement vector is replaced by a periodic 
function of the finite differences of the displacement vec- 
tor that has the same period as the lattice and is linear 
for small differences. The potential energy is obtained 
from the resulting strain energy density by summing over 
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all lattice points [8]. The equations of motion obtained 
from the crystal potential energy reduce to those of linear 
anisotropic elasticity far from the dislocation cores, where 
the differences of the displacement vector are small and 
approximately equal to their partial derivatives. At the 
dislocation cores, the singularities of the stress field are 
regularized by the lattice and gliding is allowed by the 
above mentioned periodic functions. Extensions to bec 
and fee crystals [8] or to lattices with a two-atom basis [9] 
are possible. While these models are three-dimensional, it 
is possible to simplify them for simple geometrical config- 
urations such as planar edge dislocations having parallel 
Burgers vectors. In this case, the component of the dis- 
placement vector parallel to the Burgers vector suffices to 
describe depinning and gliding of dislocations. A simpli- 
fied model retaining only this component has been used 
to analyze dislocation depinning and motion at the Peicrls 
stress [10]. The same phenomena occur in the same way 
in the more complete models having two nonzero compo- 
nents of the displacement vector [8] . 

Model. — Here we study a simple 2D periodized dis- 
crete elasticity scalar model that may describe homoge- 
neous nucleation of edge dislocations with parallel Burgers 
vectors and is amenable to detailed analysis. Consider a 
2D simple cubic lattice with lattice constant normalized 
to 1, lattice points labelled by indices i = 1, ...,N X 

and j = l,...,N y , and displacement vector (uij,0). At 
the boundary, a shear strain F is applied, so that the 
displacement Uij for j = l,N y and for j = 1,N X is 
F[j — (N y + l)/2]. F is also the dimensionlcss shear 
stress. To visualize more easily the lattice, we will some- 
times use coordinates (x,y) with x = i — (N x + l)/2 and 
y = j — (N y + l)/2 centered at the lattice center. In these 
coordinates, the boundary displacement adopts its usual 
form Fy. Uij obeys the following nondimensional equa- 
tions: 
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+A {g a {ui,j+i - Ui,j) + 9a(ui,j-i - Ui,j)l- (1) 

Here A = C44/C11 provided we consider cubic crystals 
with elastic constants Cn, C12, C44. If we select a nondi- 
mensional time scale C\\tj (pl 2 ^j) — » t, then a = 1 and 
m = Cn/ (pl 2r y 2 ), where 7 is a friction coefficient with 
units of frequency, p is the mass density and I the dimen- 
sional lattice constant. The nondimensional displacement 
vector is measured in units of I. With this choice of scales, 
we can consider the overdamped case with m = 0. On 
the other hand, if we select a nondimensional time scale 
CKS/Qp 1 / 2 ) -» t, then m = 1 and a = l^y/p/C^. With 
this second choice of scales, we can consider the conserva- 
tive case with a — 0. 

The nonlinear function g a is periodic, with period equal 
to the lattice space and g' a (0) = 1. It allows gliding of half 
columns of atoms in the x-direction, which is the direction 
of the Burgers vectors of the edge dislocations that will be 



nucleated. Gliding along other directions is not possible 
in this model: we would need a two-component displace- 
ment vector and a periodic function of finite differences 
along the x axis [8] . We have used in our simulations the 
continuous one-parameter family of periodic functions 



, x 2a f sin(ff) , -a < x < a, 

9a(x) = -\ a<x<l-a, 



(2) 



with < a < 1/2 and period 1. In the symmetric case 
a = 1/4, (1) - (2) is the interacting atomic chains model 
[11] . The parameter a controls the asymmetry of g a , which 
in turn determines the size of the dislocation core and the 
Peierls stress needed for a dislocation to start moving [8]. 
As a increases, the interval over which g' a (x) > increases 
at the expense of the interval over which the slope of g a is 
negative. Then as a increases, so does the Peierls stress, 
whereas both the core size and the mobility of defects 
decrease. Large values of a result in very narrow cores and 
large Peierls stresses 1 . The value of a can be selected so 
that the Peierls stress calculated from (1) fits experimental 
values or values calculated using molecular dynamics [8]. 

We consider first the overdamped case, m = 0, a = 1 
with A = 0.3071 (corresponding to tungsten), a = 0.2. 
In this case and for time independent shear, the potential 
energy 



1 



{ui+i,j - U -Lj) 2 + AG a (u i+ ij - Uij) 



. (3) 



with G' a (x) = g a (x) and G a (0) = is a Lyapunov func- 
tional of the gradient system (1): it satisfies V > and 

— - dV dui 'i _ ( dui -i \ 2 < n 

d± ~ ^ d Ui j dt - ' 



dt 



since duij/dt = —dV/dui_j and the shear strain F does 
not depend on time. In the unstressed crystal configu- 
ration F — 0, a given initial condition evolves exponen- 
tially fast to a stable homogeneous dislocation-free station- 
ary state which we call BR0. This stable solution of the 
discrete model is simply the undisturbed lattice without 
any dislocations. As we select larger and larger positive 
stresses, the homogeneous stationary configuration BR0 
is strained but it continues to be stable and dislocation- 
free until a critical stress F c = 0.2193 is reached. At F c , 
the maximum eigenvalue of Eq. (1) linearized about the 
stationary solution BR0 becomes zero. What happens for 
F > F c l 

Multistable configurations with dislocations. 

One way to proceed is to start from the stable station- 
ary configuration BR0 at F = 0. We then increase the 
shear strain F to a small value AF, use the configuration 



1 Note that the parameter a used in [8] corresponds to —a + 1/2 
in (2) and therefore the Peierls stress in Figure 2 of [8] decreases as 
a increases. 
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Fig. 1: Configurations of the stationary solutions (a) BR1, (b) 
BR2, at F f = 0.22. Parameter values are a = 0.2, A = 0.3071 
for a 10x10 lattice. The crosses represents the positions of 
the boundary atoms which are fixed by the shear boundary 
condition. 



BR0 for F = as initial condition, solve (1) and find the 
corresponding stable stationary configuration. Repeating 
this procedure, we follow BR0 until F c and for F > F c we 
obtain the configuration BR2 at the corresponding value 
of F depicted in Fig. 1(b). The stationary configuration 
contains four edge dislocations, two with Burgers vector 
(1,0), the other two with Burgers vectors (—1,0). These 
dislocations appear as the result of the nucleation of two 
edge dipole dislocations at y — ±1. Immediately after 
they are created, these dipoles split into their component 
dislocations that move in opposite directions until they 
reach the sample boundary. Why do dipoles split in op- 
posite moving edge dislocations? The strain needed for a 
given edge dislocation dipole to split into its component 
dislocations is about F c /10 [12], much smaller than the 
strain F c needed to nucleate a dipole. Thus the dipoles 
split into opposite moving dislocations immediately after 
being created. 

Can we obtain other configurations by doing things dif- 
ferently? The answer is yes. Suppose we want to explore 
the stable stationary configurations at a strain Ff — 0.22 
slightly larger than F c . Starting with BR0 at F = 0, we 
turn in the strain according to a linear law during a time 
t r (the ramping time) and leave F = Ff for t > t r . Then 
F(t) = c r tH{t r - t) + F f H(t - t r ), where c r = F f /t r 
and H(x) are the strain rate and the Heaviside unit step 
function, respectively. Same as in other multistable sys- 
tems [13], we obtain different final stable configurations 
depending on Ff and t r . For long ramping times t r > 87, 
we again find BR2. Fig. 2(c) and (d) show two snapshots 
of the strain component 2ei,2 = ga{ui.j+i — taken 
after F has reached its final value Ff and Uij is evolving 
towards its final stationary configuration. We observe two 
depressions of the strain ei.2 at y = ±1 indicating nucle- 
ation of two dislocation dipoles. As we have said before, 
the Peierls stress needed to split and move the component 
dislocations in a dipole is much smaller than the stress 
required for homogeneous nucleation of one dipole. Thus 
after being nucleated, the edge dislocations with oppo- 
site Burgers vectors comprising each dipole immediately 



Fig. 2: Four snapshots of the strain 2ei,2 at times (a) 1491.8 
and (b) 2643.2 for the evolution towards BR1 at ramping time 
U = 86 (cr = 0.0026), and at times (c) 2056 and (d) 3242.3 for 
the evolution towards BR2 at ramping time t r — 1000 (c r = 
2.2 x 10 -4 ). Ff = 0.22 and other parameter values as in Fig. 1. 



move in opposite directions towards the lattice boundaries. 
The final configuration is BR2. For shorter ramping times 
(82 < t r < 87), Figs. 2(a) and (b) show that only one dislo- 
cation dipole is nucleated, splits into two edge dislocations 
with opposite Burgers vectors that then move towards the 
boundaries in opposite directions. The final configuration 
is BR1 as in Fig. 1(a). For t r < 81, a final configuration 
BR3 (similar to that in Fig 1(b) but not explicitly shown) 
is reached after two dipoles are nucleated at the upper and 
lower boundaries and their component dislocations move 
to the left and right boundaries in opposite directions. 

Nucleation time. — Fig. 3 shows the evolution of 
the potential energy when the same final strain Ft is 
reached at three different strain rates. We observe that 
for very short ramping times, e.g., t r = 10, the energy 
reaches a peak (with strain energy larger than that of the 
homogeneous branch Vbro = 0.6335) at t r and then re- 
laxes toward its final value Vrr3 = 0.6433 correspond- 
ing to configuration BR3. For very long ramping times, 
V(t) follows BR0 adiabatically beyond t r , at the plateau 

V = Vbro that lasts from t r until 1t T . At this later time, 
the energy drops abruptly to its final value correspond- 
ing to configuration BR2 in Fig. 1(b) with strain energy 

V = Vbr2 = 0.5417. The abrupt energy drop marks 
the nucleation of the two dipoles. Similarly, the precip- 
itous energy drop in the evolution towards configuration 
BR1 of Fig. 1(a) with energy Vrri = 0.5143 corresponds 
to the nucleation of one dipole after a long plateau with 
strain energy Vbro ends abruptly at time 17 t r . Note that 
the evolution towards BR1 starts with a small spike at 
t = t r (for an intermediate ramping time of 86), it contin- 
ues with energy V — 0.6435 at a short plateau for times 
1 < tjt T < 2.5, there is a gradual energy decay to the 
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Fig. 3: Energy V(t) for ramping times 10 (BR3), 86 (BR1) and 
1000 (BR2). Parameter values are the same as in Fig. 1. 



long plateau at Vbro = 0.6335 that lasts from t = 3.5 t r 
till t = 17t r , and then the strain energy drops to its final 
value Vbri = 0.5143. 

Bifurcation diagram. — More precise information 
about possible stable configurations can be obtained by 
examining the bifurcation diagram of the I 2 norm of the 

displacement, ||u|| = \ Yl, u 'i i (the sum excludes points at 



the boundaries), versus the strain F 2 The complete bifur- 
cation diagram has been calculated using the AUTO pro- 
gram of numerical continuation of solutions [14], and it is 
rather complex: there are many bifurcation points issuing 
from different stationary solution branches, most of which 
are unstable. If we depict all possible solution branches, 
the resulting bifurcation diagram is rather messy. Thus 
we have chosen to depict only important solution branches 
which are stable in certain strain intervals. 

In Fig. 4(a) we show the only two primary branches that 
bifurcate from BR0 at F = F c = 0.2193. The inset shows 
that these branches appear as a subcritical pitchfork bifur- 
cation at F — F c . Both start being unstable for F close to 
F c but become stable after limit points (BR1 exactly after 
the limit point, BR2 becomes stable after a secondary bi- 
furcation point with F > F12), giving rise to intervals were 
several stationary solutions are simultaneously stable. As 
explained before, these configurations can be selected by 
turning the final shear strain at different rates. 

The branches BR1 and BR2 contain a number of 
secondary bifurcations and limit points, as depicted in 
Fig. 4(b). These branches fold over themselves in segments 



2 We could have depicted the potential energy V versus strain 
as a bifurcation diagram. However, the energy of the bifurcating 
branches BR1 and BR2 near the main subcritical bifurcation point 
is so close to the energy of the branch BR0 that visualizing this 
bifurcation is very hard in an energy-strain diagram. Thus we have 
preferred to use the P norm. 
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Fig. 4: (a) Bifurcation diagram showing only the primary sta- 
tionary branches issuing from the homogeneous solution BR0. 
At F c , branches BRI and BR2 appear as a subcritical pitch- 
fork bifurcation from BR0 (see the insets), (b) Larger view of 
the bifurcation diagram showing branches that are prolonga- 
tions of BRI and BR2 after a number of limit and bifurcation 
points. Solid lines correspond to stable solutions, dashed lines 
to unstable solutions, limit points are marked as triangles and 
bifurcation points as squares. Parameter values are the same 
as in Fig. 1. 



delimited by additional limit points and display other lin- 
early stable parts. The configurations thereof contain ad- 
ditional dislocations arising from dipole nucleation. For 
example, BRI has another stable configuration at Ff aris- 
ing from nucleation of dipoles at y = ±2, whereas BR2 
has stable configurations arising from dipole nucleation at 
y = 0, ±3 and at y — ±3, respectively. While these con- 
figurations are linearly stable, we have not been able to 
reach them by ramping to Ff — 0.22. 

It is interesting to note that the strains 2ei.2 correspond- 
ing to the unstable parts of BRI and BR2 resemble the 
snapshots corresponding to dipole nucleation in Fig. 2. If 
we follow the unstable part of BRI backwards from the 
limit point at Fn ~ 0.11 to the critical strain F c , we ob- 
serve that its strain 2ei ; 2 has a depression at x = cor- 
responding to dipole nucleation for Fu < F < F c , but 
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this depression becomes less and less observable as we ap- 
proach F c . Similarly, as we follow the unstable part of 
BR2 from its limit point Fi 2 ~ 0.17 to F c , we observe that 
2ei j2 first exhibits two symmetric depressions at x = cor- 
responding to nucleation of two dipoles. Then these de- 
pressions diminish until the configuration of the unstable 
part of BR2 becomes very similar to that of BR1 as F ap- 
proaches F c . This is as it should be because BR1 and BR2 
merge at F c in a subcritical pitchfork bifurcation: Near F c , 
BR1 and BR2 differ from BRO by cx ±y/F c - Ftpij 
is the eigenvector corresponding to the zero eigenvalue at 
F = F c ). 

The components of the eigenvector ipij have oppo- 
site signs for alternate rows of their corresponding lattice 
sites 3 , i>i,ji>i,j+i < 0, while they keep the same sign along 
the same column in the lattice, Vi,jV'i,fe > 0. The differ- 
ences \tpij-\-i — ipi,j\ are found to be largest at the center 
of the lattice. We have observed that regions where dislo- 
cations may nucleate are located between rows j and j + 1 
for which ipij+i — ipij is maximum (y = 0, ±2 for branch 
BR1) or minimum (y = ±1,±3 for branch BR2). 

Influence of lattice size on bifurcations. In all 

our figures, we have presented numerical solutions corre- 
sponding to a 10x10 lattice, a = 0.2 and A = 0.3071. The 
secondary solution branches (not depicted here) change 
substantially with N x , N y , A and a. Bifurcation and 
limit points also might appear or disappear from the di- 
agram. However, BR1 and BR2 persist and the value of 
F c and the type of the primary bifurcation do not change 
if we increase the computational domain. For example, 
F c = 0.2053 for a 20x20 lattice which is closer to F c = a 
than the critical strain for a smaller lattice: apparently 
F c — > a as the lattice size increases. In some cases (small a 
and small lattices, such as a = 0.1 in 6x6 and 8x8 lattices), 
the pitchfork bifurcation at F c is supercritical. However, 
the bifurcation becomes again subcritical for larger val- 
ues of a and for larger lattices (a = 0.1, in 10x10 and 
14x14 lattices). Finding the branch BR3 with AUTO has 
been quite elusive and, in fact, we did not find it for the 
10x10 lattice with the parameter values of Fig. 1 using 
AUTO, whereas ramping produced BR3 in a straightfor- 
ward manner. For a 6x6 lattice with a = 0.4, we find that 
BR3 appears as one of the branches issuing from BRO as 
a supercritical pitchfork bifurcation at F = F 3 > F c . This 
branch has a limit point at a larger strain and it continues 
for decreasing values of F. One of the stretches of BR3 is 
linearly stable at a range of F overlapping those of BR1 
and BR2 [12]. 

Effects of inertia. — The bifurcation diagram corre- 
sponding to stationary solutions of (1) with inertia (m ^ 
0) is still the same as presented here. However, the stabil- 
ity character of the solutions changes. In the conservative 

3 Note that lattice sites having the same coordinate j (y axis) 
form a row and sites having the same coordinate i (x axis) form a 
column. This is different from the usual convention denoting rows 
and columns of a matrix A; j. 



case m = 1, a = 0, stable solutions are no longer asymp- 
totically stable. Linearizing (1) about a stable solution, we 
find a problem with purely imaginary eigenvalues. There- 
fore these solutions are centers: small disturbances about 
them give rise to small permanent oscillations about them. 
The linearized problem about unstable solutions has pairs 
of positive and negative eigenvalues and therefore these 
solutions are saddle-centers in general. 

Concluding remarks. — What have we learned 
about homogeneous dislocation of nucleations by shear- 
ing a dislocation-free state? Clearly the critical strain 
F c marks the instability of the dislocation-free solution 
branch BRO. F c is characterized as the shear strain at 
which the largest eigenvalue of the linear eigenvalue prob- 
lem about BRO becomes zero. The components of the cor- 
responding eigenvector indicate possible nucleation sites 
that are realized by different stationary solution branches 
in the bifurcation diagram. The fact that the pitchfork bi- 
furcation at F c is subcritical implies that dislocation nucle- 
ation can occur at subcritical shear strain values at which 
the solution branches BR1 and BR2 become stable. Thus 
it is important to determine the ranges of F at which some 
of the branches BRO, BR1 and BR2 are simultaneously 
stable. This cannot be done by simple linear stability cal- 
culations: instead numerical continuation algorithms such 
as AUTO have to be employed. For overdamped dynam- 
ics, different stable stationary configurations can be se- 
lected by the strain rate at which the final strain Ff is 
reached. An abrupt drop in energy marks the nucleation 
time at which one or two dipoles are nucleated. Since the 
Peierls stress is lower than the critical stress for homo- 
geneous dipole nucleation, the dipoles immediately split 
into edge dislocations with opposite Burgers vectors mov- 
ing in opposite directions. This bifurcation picture seems 
to describe larger lattices and it captures the stationary 
solutions even if inertia is added. 

Experimental studies often use ad hoc criteria for nu- 
cleation of dislocations such as the critical resolved shear 
stress (CRSS) [1]. To use this criterion, the critical stress 
for nucleation has to be related to the applied force by 
other means, such as the Hertz contact theory in nanoin- 
dentation experiments [1,15]. Moreover, the critical stress 
itself has to be calibrated independently and it cannot 
be a fixed value. Instead, the ideal shear stress for nu- 
cleation may depend strongly on the other stress com- 
ponents, not just on the shear stress component acting on 
the plane, as shown by density functional theory [16]. Ear- 
lier continuum mechanics studies suggest that nucleation 
of dislocations is related to the loss of strict convexity in 
the energy and stress concentration [17]. Recent studies 
calculate the elastic constants and internal stresses from 
atomistic calculations or from finite element calculations 
and the Cauchy-Born hypothesis to figure out atom mo- 
tion [18]. Then they minimize a certain scalar functional 
of elastic constants and internal stresses at each point of 
the solid. Nucleation occurs at those points at which the 



p-5 



I. Plans et al. 



resulting scalar functional first vanishes [19]. There is a 
widespread feeling in all these studies that homogeneous 
nucleation of dislocations is related to some bifurcation oc- 
curring once the instability starts but they do not report 
any precise analysis and calculation of this bifurcation, in 
contrast to our work. 

Dislocation depinning and motion and dislocation inter- 
action occur in the same way in the simple scalar model 
(1) [10] and in more complete planar discrete elasticity 
models with two components of the displacement vec- 
tor [8] . Thus we expect that our bifurcation description of 
homogeneous nucleation and motion of dislocation dipoles 
also applies to these planar models. Studies of nucleation 
of dislocations in more complete two and three dimen- 
sional models are postponed to future publications. 
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